/**
 * ANSWER: 19823.542204
 */

#include <iostream>
#include <vector>
#include <numeric>
#include <algorithm>

// Solve the linear system:
// x_i = x_i
// x_k = a[k]*x[k-1] + (1-a[k])*x[k+1]
// x_j = j
static void solve_subsystem(int n, int i, int j, const std::vector<double> &a, std::vector<double> &x)
{
	static std::vector<double> p, q;
	p.resize(n);
	q.resize(n);

	p[i] = 0;
	q[i] = x[i];
	for (int k = i + 1; k < j; k++)
	{
		p[k] = (1-a[k])/(1-a[k]*p[k-1]);
		q[k] = a[k]/(1-a[k]*p[k-1])*q[k-1];
	}

	for (int k = j - 1; k > i; k--)
	{
		x[k] = p[k]*x[k+1] + q[k];
	}
}

#if 0
void solve_problem_339()
{
	const int N = 100; // total number of sheep

	int cc = 0;
	std::vector<double> x, a, last_x;
	std::vector<bool> bumped;
	x.push_back(0);
	
	for (int n = 1; n <= 2*N; n++)
	{
		last_x.assign(x.cbegin(), x.cend());

		x.resize(n+1);
		a.resize(n);
		bumped.resize(n);
		x[n] = n;
		cc = 0;

		// Initialize weights a[k] and bumping vector
		for (int k = 1; k < n; k++)
		{
			a[k] = (double)(n-k) / (double)n;
			bumped[k] = false;
		}

		// For each unbumped x[k], check whether x[k] >= a[k]*x[k-1] + (1-a[k])*x[k+1].
		// If not, mark x[k] as "need to bump".
		bool do_loop = true;
		while (do_loop)
		{
			do_loop = false;
			int last_unbumped = n;
			for (int k = n - 1; k >= 1; k--)
			{
				if (!bumped[k])
				{
					double t = a[k]*x[k-1] + (1-a[k])*x[k+1];
					if (x[k] < t)
					{
						do_loop = true;
						bumped[k] = true;
						x[k] = t;
					}
					else
					{
						if (last_unbumped > k + 1)
						{
							solve_subsystem(n, k, last_unbumped, a, x);
							cc++;
						}
						last_unbumped = k;
					}
				}
			}
			if (last_unbumped > 1)
			{
				solve_subsystem(n, 0, last_unbumped, a, x);
				cc++;
			}
		}
		
		//if (n % 200 == 0)
		//{
		//	std::cout.precision(6);
		//	std::cout << "E(" << n/2 << ") = " << std::fixed << 0.5*(x[n/2-1]+x[n/2+1]) << std::endl;
		//}
		std::cout << "Free variable = " << (n-1);
		std::cout << ", bumped = " << std::count(bumped.begin(), bumped.end(), true) << std::endl;

		if (n == 2*N)
		{
			std::cout << "Total = " << n << std::endl;
			std::cout << "Bumped = " << std::count(bumped.begin(), bumped.end(), true) << std::endl;

			// try solve it at once??
			std::vector<double> y(n+1);
			y[0] = 0;
			y[n] = n;
			solve_subsystem(n, 0, n, a, y);

			// how many over last?
			int many = 0;
			for (int k = 1; k < n; k++)
			{
				if (y[k] > last_x[k])
					many++;
			}
			std::cout << "At once? bumped = " << many << std::endl;
		}
	}

	//for (int k = 0; k <= 2*N; k++)
	//{
	//	std::cout << "f(" << k << ", " << (2*N-k) << ") = " << x[k] << std::endl;
	//}
	std::cout.precision(6);
	std::cout << "E(" << N << ") = " << std::fixed << 0.5*(x[N-1]+x[N+1]) << std::endl;
	std::cout << "Called solver " << cc << " times." << std::endl;
}
#else

void solve_problem_339()
{
	const int N = 10000; // number of sheep in each flock

	std::vector<double> x, a;
	x.push_back(0);
	
	for (int n = 1; n <= 2*N; n++)
	{
		x.resize(n+1);
		a.resize(n);
		x[n] = n;

		// Initialize weights a[k] and bumping vector
		for (int k = 1; k < n; k++)
		{
			a[k] = (double)(n-k) / (double)n;
		}

		// We have the following unproved rule: 
		// - the mid point is unchanged; 
		// - all points to the left are unchanged;
		// - all points to the right need to be solved.
		solve_subsystem(n, n/2, n, a, x);

		//if (n % 200 == 0)
		//{
		//	std::cout.precision(6);
		//	std::cout << "E(" << n/2 << ") = " << std::fixed << 0.5*(x[n/2-1]+x[n/2+1]) << std::endl;
		//}
		if (0)
		{
			std::cout << "[" << n << "] ";
			for (int k = 0; k <= n; k++)
			{
				std::cout << x[k] << " ";
			}
			std::cout << std::endl;
			if (n == 10)
				return;
		}
	}

	std::cout.precision(6);
	std::cout << std::fixed << 0.5*(x[N-1]+x[N+1]) << std::endl;
}

#endif
